function [GIRFxhr,GIRFy,GIRFx] = getGIRFs(model,theta2,x0,z0,numSim,IRFlength,shockSize)

% Unfold the factor dynamics
nx = size(model.gx,2);
[h0_1,h0_2,hx_1,hx_2,sigma,gama0,gamaz,gamax,sigmaz]   = unfoldTheta2_threshold_macro(theta2,nx);

GIRFxhr   = zeros(length(model.A),model.nx,IRFlength);
GIRFy     = zeros(model.ny,model.nx,IRFlength);
GIRFx     = zeros(model.nx,model.nx,IRFlength);

rng(1,'twister')
epsilonx                     = zeros(numSim,IRFlength,nx);
epsilonx(1:numSim/2,:,:)     = randn(numSim/2,IRFlength,nx);
epsilonx(numSim/2+1:end,:,:) = -epsilonx(1:numSim/2,:,:);

epsilonz                      = zeros(numSim,IRFlength,1);
epsilonz(1:numSim/2,:,:)      = randn(numSim/2,IRFlength,1);
epsilonz(numSim/2+1:end,:,:)  = -epsilonz(1:numSim/2,:,:)   ;
for i=1:nx
    i
    for s=1:numSim
        xs_cu = x0;
        zs_cu = z0;
        x_cu  = x0;
        z_cu  = z0;
        
        for k=1:IRFlength
            % Path with shock
            epsilonxUsed      = squeeze(epsilonx(s,k,:));
            if k == 1
                epsilonxUsed(i,1) = shockSize;
            end
            if z_cu >= 0
                xs_cup = h0_1 + hx_1*xs_cu + sigma*epsilonxUsed;
            else
                xs_cup = h0_2 + hx_2*xs_cu + sigma*epsilonxUsed;
            end
            zs_cup = gama0 + gamaz*zs_cu + gamax'*xs_cu + sigmaz*epsilonz(s,k,1);
            
            % Path without shocks
            epsilonxUsed      = squeeze(epsilonx(s,k,:));
            if z_cu >= 0
                x_cup = h0_1 + hx_1*x_cu + sigma*epsilonxUsed;
            else
                x_cup = h0_2 + hx_2*x_cu + sigma*epsilonxUsed;
            end
            z_cup = gama0 + gamaz*z_cu + gamax'*x_cu + sigmaz*epsilonz(s,k,1);
            
            % The observables with shock
            yields_s = model.g0 + model.gx*xs_cup + model.gxx*kron(xs_cup,xs_cup);
            xhpr_s = getExpectedExcessReturns_h1(model,theta2,xs_cup,zs_cup)*12;
            
            % The observables without shock
            yields = model.g0 + model.gx*x_cup + model.gxx*kron(x_cup,x_cup);
            xhpr = getExpectedExcessReturns_h1(model,theta2,x_cup,z_cup)*12;
            
            % The GIRFs
            GIRFxhr(:,i,k) = GIRFxhr(:,i,k) + xhpr_s-xhpr;
            GIRFy(:,i,k)   = GIRFy(:,i,k)   + yields_s - yields;
            GIRFx(:,i,k)   = GIRFx(:,i,k)   + xs_cup - x_cup;
            
            % Updating the states
            xs_cu = xs_cup;
            zs_cu = zs_cup;
            x_cu  = x_cup;
            z_cu  = z_cup;
            
        end
    end
end
GIRFxhr = GIRFxhr/numSim;
GIRFy   = GIRFy/numSim;
GIRFx   = GIRFx/numSim;
% For GIRFxhr, reporting same maturities as for GIRFy
GIRFxhr = GIRFxhr(model.matSelect,:,:);

